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Abstract 

We present a novel formulation for biochemical reaction networks in the context 
of protein signal transduction. The model consists of input-output transfer functions, 
which are derived from differential equations, using stable equilibria. We select a set 
of 'source' species, which are interpreted as input signals. Signals are transmitted to 
all other species in the system (the 'target' species) with a specific delay and with a 
specific transmission strength. The delay is computed as the maximal reaction time 
until a stable equilibrium for the target species is reached, in the context of all other 
reactions in the system. The transmission strength is the concentration change of the 
target species. The computed input-output transfer functions can be stored in a 
matrix, fitted with parameters, and even recalled to build dynamical models on the 
basis of state changes. By separating the temporal and the magnitudinal domain we 
can greatly simplify the computational model, circumventing typical problems of 
complex dynamical systems. The transfer function transformation of biochemical 
reaction systems can be applied to mass-action kinetic models of signal transduction. 
The paper shows that this approach yields significant novel insights while remaining a 
fully testable and executable dynamical model for signal transduction. In particular 
we can deconstruct the complex system into local transfer functions between 
individual species. As an example, we examine modularity and signal integration 
using a published model of striatal neural plasticity. The modularizations that 
emerge correspond to a known biological distinction between calcium-dependent and 
cAMP-dependent pathways. Remarkably, we found that overall interconnectedness 
depends on the magnitude of inputs, with higher connectivity at low input 
concentrations and significant modularization at moderate to high input 
concentrations. This general result, which directly follows from the properties of 
individual transfer functions, contradicts notions of ubiquitous complexity by showing 
input-dependent signal transmission inactivation. 



Introduction 

Biochemical reaction systems are usually conceptualized as dynamical systems - systems 
that evolve in continuous time and may or may not receive additional input to the system. 
Mathematically, this can be expressed by sets of ordinary differential equations (ODE), such 
that rates of concentration changes correspond to mass-action kinetic parameters [T8l IT9] . 
In this paper we use existing mass-action dynamical systems to propose an alternate or 
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additional framework for modeling and interpretation of biochemical reaction systems. We 
provide an algebraic analysis of biochemical reaction systems as a matrix of concentrations for 
all species, given certain input concentrations. These concentrations correspond to steady- 
state amounts which are reached after a delay time, and the delay times can be measured 
by the system as well. 

We use an arbitrary published model p] as an example for a ODE dynamical model of 
biochemical reactions. The model simulates intracellular signal transduction from receptor 
binding to molecular targets in different cellular compartments, as an important component 
in the long-term regulation of protein expression implied in neural and synaptic plasticity. 
In striatal neurons, both a calcium-dependent pathway and a cAMP-dependent pathway are 
activated during the initiation of neural plasticity by NMDA/AMPA receptors and neuro- 
modulator receptors such as dopamine Dl receptors [2}|3lll]. Their effects and the integration 
of signaling on common targets such as kinases and phosphatases have been the subject of 
a number of computational models [5J El El El E] . In particular, the role of the DARPP- 
32 protein in striatal neurons in determining the outcome of membrane signaling has been 
modeled by different groups, based on a common set of experimental data [T0| fTT| [T], cf. 
[T2l [TU [151 HEl [17]. Many similar models [23] have been developed in the last 10-12 years in 
different areas of biology. Models with dozens or more of species have up to a 100 or more 
equations and are consequently complex and difficult to understand as continuous dynamical 
systems |13j . A transformation into a matrix-based formulation of input-output functions, 
even at the cost of a loss of fast dynamical modeling, promises considerable gain of insight 
and access to a different set of mathematical tools. Simple mass-action kinetic models may 
be criticized for disregarding the real complexity of spatio-temporal molecular interactions. 
Some alternatives use spatial grids and stochastic versions of biochemical reactions to cap- 
ture this complexity [201 121] • However, certain variations, such as compartmental modeling 
with diffusion, altered kinetics for anchored proteins, or employing molecular kinetics as the 
basis for binding constants may be employed within the mass-action kinetic framework to 
achieve better correspondence with the biological reality. These variations can be directly 
transferred to the proposed model as well. 

In our approach, we identify input nodes, and then pre-compute the outcomes for all in- 
ternal species (target species) in response to biological meaningful ranges and combinations 
of inputs. This allows to analyze a biochemical reaction system under all possible input 
conditions. The analysis can be done for arbitrary ODE models [23], provided minimal re- 
quirements on conservation properties are realized (cf. section "Methods", [21] )■ The results 
are stored as vectors or matrices ('systemic protein signaling functions' (psfs)) and can be 
fitted with functional parameters. It is an important aspect of the model that computa- 
tions are done systemically. In section "Elementary Biochemical Reactions", source-target 
interactions are first analyzed in isolation ('elementary psfs'). They all constitute hyperbolic 
saturation functions, therefore rate parameters can be uniformly translated into functional 
parameters for signal transmission strength. But in a systemic context, source-target inter- 
actions change because of additional influences on the species from other equations in the 
system (section "Systemic PSF Analysis"). Therefore a fitted systemic psf from A to B is 
different from the elementary psf. The pre-computed, systemic psfs may be used to cre- 
ate state-change simulation models, i.e. discrete-time models, which can be compared with 
continuous ODE models (section "Systemic Delay and State-change Dynamics"). What is 
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significant and novel about our analysis is that we can extract systemic transfer functions 
from the complex system, and thereby dissect the system into parts. We can analyze the 
transmission properties of individual species, compare their minimum and maximum values, 
and the functional shape of their transmission strength. Specifically, we can show under 
which circumstances a link is functional, i.e. actually transmits information (section "Com- 
puting Input-dependent Modularity"). 

The analysis has a number of restrictions. An important restriction is that our model 
does not allow for analysis of fast interactions below the resolution of settling into steady- 
state. The requirement of conservation of mass guarantees that for each input concentrations 
will eventually settle to some equilibrium value, but due to the prevalence of feedback in- 
teractions, they may still produce transients or dampened oscillations. This means that 
fast fluctuations of input will not be adequately simulated using pre-computed psf functions 
alone. It is then necessary to refer to the underlying ODE model. The model is most suitable 
for studying disease states, pharmacological interventions, genetic manipulations, miRNA in- 
terference, or any system conditions which fundamentally alter the presence or concentration 
of molecular species. These conditions may then be tested either in steady-state or with a 
sequence of sufficiently slow input constellations. A second restriction is that the model in- 
herits parameter uncertainty from mass-action kinetic models. These parameters are derived 
from experimental measurements, but typically with a high degree of uncertainty [22]. Our 
analysis offers a clear distinction between elementary and systemic functional parameters 
and explains why experimental measurements are so highly dependent on the systemic con- 
text. In this paper, we have treated elementary parameters as given by the underlying model 
[I]. In principle, systemic functional interactions can be measured experimentally, and in 
the model each interaction can be adjusted separately. This may offer a novel theoretical 
approach towards finding adequate elementary rate parameters. 

Materials and Methods 
System Definition 

A biochemical reaction system formulation for signal transduction contains two different 
types of reactions: 

1. complex formation 
[A] + [B] o [AB] 

2. enzymatic reactions 

[A] + [E] o [AE] -> [A*] + E. 

The system has concentrations for species A, B and E (A*, AE and AB can be calculated), 
a set of kinetic rate parameters k on and k /f for the forward and backward binding reactions 
in complex formation, and k cat for the rate of enzymatic production. The equational struc- 
ture, the kinetic parameters and the initial concentrations of the model pp are reproduced in 
Tabled]- Table HJ with slight modifications: Equation 40 was added to 'close the loop' from 
AMP to ATP and thus provide for conservation of all molecules in the system. Conservation 
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Table 1. Reactions in the cAMP pathway 



of mass is necessary in order for all species to reach equilibrium. It means that for any for- 
ward reaction there needs to be a reverse reaction, such that any species receives both input 
and output ('weakly reversible system', cf. [21] )• This implies that pure loss reactions, like 
endocytosis or diffusion across the cell membrane (secretion) cannot adequately be modeled 
with this system, unless balancing reactions are added which make up for the loss. E.g., 
for endocytosis of a ligand-bound receptor, both the ligand and the receptor, possibly in- 
dependently, have to be recycled, which means that bridging equations for receptor-ligand 
dissolution in the endocytosed state, and input rates for receptors and ligands have to be 
added. Secondly, the species PP1 and its interactions (4 equations) were left out, since they 
contain a complex which dissolves into its 3 components in one step: this would require a, 
fairly trivial, addition to the current psf system implementation. The system can be depicted 
as a bipartite graph with nodes for species and nodes for reactions. 

PSF Analysis 

In order to set up source-target functions, we need to select input nodes from the available 
species nodes. In this example, we used Da (dopamine as ligand for the Dl receptor) and 
Ca (extracellular calcium that diffuses through ion channels in the membrane). We use 
input concentrations over a specified range (e.g., between 60nM and 5/iM for Da), sample 
over the range with e. g. n=20 steps, and use the differential equation implementation of 
the system to calculate the output values for all species for each sampling step. Because of 
the conservation of molecules, all species reach steady-state after a sufficient period of time. 
We define steady-state pragmatically by relative change of less than 2% over 100s. We also 
use the established terminology of EC10, EC50, EC90 etc. to indicate 10, 50 or 90 % of 
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Table 3. DARPP-32 reactions 



steady-state concentration value. Additionally, we calculate the delay in reaching steady- 
state. We store input-target concentration mappings in a vector (single-input system), or a 
matrix (multiple-input system). We fit the vectors with hyperbolic or linear functions, using 
standard techniques in Matlab (fminsearch, [25]). In this way we derive parameters which 
can be analyzed and used instead of the explicit vectors. (In this paper, the fitting is only 
done for single-input systems, multiple-input systems require different techniques.) 

All information on source-target transfer functions for the complete, complex signaling 
system ('systemic psf ') can be stored in a static data structure. For each species, it contains 
its concentration range, and for each reaction, it contains the parameters of the functional 
fit. We gain the possibility to regard any species as source and any other species as target 
(they may be coupled by an arbitrary number of reactions) and obtain a systemic psf as the 
transfer function between them. This representation allows to analyze the complex signaling 
system by its parts, i.e. as a set of matrices or vectors, which is the main achievement 
relative to the ODE dynamical system. In addition, dynamical simulation with appropriate 
update times may be realized by the psf representation alone, i.e. the psf simulation is not 
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Initial Concentrations 



in itself atemporal, but only discrete and fairly slow. 

We visualize the data structure as a bipartite graph, and label it with the calculated 
numeric values. Each species node is labeled with its attainable concentration range given the 
input range. For complex formation reactions, we show both [A] — > [AB] and [B] — > [AB]. 
For enzymatic reactions we show [E] as the source and [A*] as the target ([E] — > [A*]). The 
result is a labeled bipartite graph, called a 'weighted dynamic graph'. 

Results 

Elementary Biochemical Reactions 

We want to represent a biochemical reaction by a time-independent signal transfer function, 
such that y = f(x) for two species x, y. We do this by designating a source species x and 
then calculating the steady-state value for another species, the target species y, for any value 
of x, given the differential equations for the biochemical reaction. For complex formation 



where the total concentrations for [A] and [B] and kinetic rate parameters kon, koff (with 
Kd = ) are given, the differential equations are: 



We may now calculate the concentration values f(x) = y for a target species [AB] given 
a range of input values for x, e.g. the source species [A]. (Fig. [Q. 

In this way we separate the calculation of the signal response magnitude, i. e. the steady- 
state concentration, from the calculation of the time until a steady-state value is reached, 
the delay. For different x, f(x) will be reached after a variable delay (Fig. [2]). 

With some modification, the same transformation applies to enzymatic reactions. The 
kinetic rate parameters are Kd = ^j^L (f or ^ anc i k cat (f or 
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exists. (For instance, 

[A*] + [E 2 \ <-> [A*£ 2 ] [A] + [£ 2 ] 
is a reaction that reverses [A*]. ) The differential equations, with k ca t2 for [A*] — > [A], are: 

dxdt(A) = fc o// [ABi] - £w[-4pi] + fc C at 2 [^*] 

dxdt{E x ) = (k off + fccat)^] - fc on [A] * [Ex] 

eted*(A*) = fccattABi] - k cat2 [A*\ 

dxdt{AEx) = k m [A][Eh] - (k off + k cat ) [AE\] 

Given concentrations for Ei, A and kinetic rate parameters Kd, k cat and k cat 2, we may 
now derive a function with x as the source species [Ei] and y as the target species [A*] 

(Fig. ED. 

In both cases the resulting curve can be fitted by a saturating hyperbolic function. 

r/\ / Vmax ymin \ 

U = j{ x ) = Umax ~ { ^ — ( — ) n 

Here y m in, the baseline concentration, is usually set to 0. 

If we choose [A] as the target of [El], we get a negative slope psf. 

r / \ ( Vmin Vmax \ 

V = J (X) = Vmin ~ ( 



! + (§)" 



We call this function the elementary protein signaling function or elementary psf. 
This function is somewhat related to a Hill equation [221 EZ]- A Hill equation is a function 
fitted to an experimental measurement to derive a dose-response relationship, comparable 
to the psf. The Hill equation allows to calculate a fractional concentration 9 for the target 
(e.g. a receptor-ligand complex) from the source concentration [L], given Kd, and fitting a 
parameter n for the steepness of the curve. 

e- W 



Kd + [L] r 



The concentration of the other compound of the complex is not used (assumed large), 
and the absolute magnitude of the target is not calculated. An equivalent for enzymatic 
reactions is not defined. The parameter n allows to measure the effect of competing binding 
reactions (n=l if none are present), which in our terminology translates into a systemic psf 
with multiple binding partners for a single target compound. Systemic psfs are a more general 
concept than Hill equations, but they relate to the same type of data, namely dose-response 
functions in steady-state. 

We have seen that signal transmission strength is uniformly characterized by saturating 
hyperbolic functions. This means that it is highest for low x and diminishes as x increases 
(Fig.HJ). For instance, in Fig. HI a 100% signal increase leads to 100%, 18% or only 6% increase 
in the target depending on the source concentration. For enzymatic reactions, absolute 
concentration changes have different effects for sources and targets of a signaling interaction. 
Signal transmission strength depends on the absolute concentration of the source, the target 
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concentration is irrelevant. This is an important observation, since protein signaling systems 
are subject to long-term regulation of concentrations. In the context of disease states or other 
sources of protein expression up-/downregulation, independence of transmission strength 
from target concentration may be an important conservative property. 

Signal transmission is strongest if a source species is expressed at a low concentration. We 
need to bear in mind however, that reaction velocity operates inversely to signal transmission 
strength: a low source concentration means a slow reaction (Fig. [2]). A functioning signaling 
system would therefore have to use an intermediate range to maximize signal transmission 
within time constraints. Our analysis opens a new way of analysis for a signaling system: 
Optimization techniques could find a best source range for both time and signal transmission 
constraints. 

Systemic PSF Analysis 

A source-target psf can be derived for any pair of species in a complex biochemical reaction 
system. For a complex system, or set of equations, we define a set of input nodes, and 
compute the output values for each possible input configuration. The analysis gives us the 
output concentration range (notwithstanding transients in a dynamic context, s. below) 
for each species, as well as a (fitted) function, or matrix of input-output correspondences. 
A biochemical reaction will produce a different psf, when it is elementary or when it is 
embedded in a context, where the participants of the elementary reaction also participate in 
other reactions. This is true for both protein complex formation and enzymatic reactions. 
We therefore call source-target functions 'systemic psfs', when they are derived from the 
context of a specific signaling system. 

We have provided this analysis for the example system. We show the concentration 
ranges and the signal transmission functions for the whole system [I] as a weighted dynamic 
graph for Da as a single input (Fig. [5]). We label each species node with its concentration 
range, determine source and target species nodes for each reaction node, and provide fits for 
the systemic psf, the transfer function that characterizes each reaction. We see from (Fig. \5§ 
that a number of systemic psfs can be fitted well with a linear function (y = mx + 6), 
showing that systemic psfs sometimes consist only of a short section from a full mapping of 
concentration values. Also, many species have only small concentration ranges, which means 
they don't have much response to Da input. 

It is an obvious advantage of the psf analysis that we are able to dissect the complex 
system and extract local properties, such as concentration ranges of individual species, and 
transfer functions for individual reactions under input stimulation. This allows to critically 
analyse a model, compare these properties with biological data, and adjust or improve the 
model in a detailed manner. 

In Fig. [HI the concentration ranges for some target species are given. We see, for in- 
stance, that among DARPP-32 phosphorylation variants, pThr75 is always more abundant 
that pThr34, by an order of magnitude. This is an example of a high-level property, which 
could be related to biological data. As another example, we notice that the active receptor 
conformation (Da-DIR) remains below 160nM even under stimulation with \\iM Da and 
more. With a D1R total concentration of 500nM, we could adjust the ligand binding coeffi- 
cient to produce more or less active receptors. Finally, the analysis shows a very low maximal 
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PKAc level (12nM) in spite of a total PKA concentration of 1.2/iM. In the original model 
PQ, blind parameter adjustment has probably generated a very low level of PKAc in order 
to achieve high signal transmission for phosphorylation of the target species pThr34, which 
is experimentally required, but which could be achieved in other ways (e.g. PP2B) as well. 

With our analysis, properties of individual species become apparent, and they can be com- 
pared to biological data, tested and adjusted on a localized basis. Even more interestingly, 
we could look for principles of 'rational system design', for instance question the transmission 
of a seven-fold increase of cAMP in the /iM range to a maximal three-fold increase in the 
10's of nM for PKAc, and analyze given biological systems from this perspective. 

In addition to the concentration ranges, we also have access to the functional mapping 
between species in the model. The systemic psfs, like the elementary psfs, are stored as 
vectors, which are matched by functional parameters. The advantage of the psf analysis is 
that we can probe a complex system on a single reaction level because the influence of the 
cellular context is encoded in the systemic psf. Thus we can compare the elementary psf with 
its transformation as a systemic psf for individual reactions. Fig. [7A shows elementary and 
systemic psfs for G-protein activation of AC5 and the calcium-activated complex AC5Ca. We 
see that the systemic psfs are somewhat deflected, compared to the elementary psf, which is 
what we expect from the parallel activation of AC5Ca and AC5 by the same species. We may 
specify a desired psf using only functional parameters, and adjust elementary parameters to 
match the psf (Fig. OB). A local change to the Kd binding coefficient between AC5Ca and 
GoaGTP allows a change in the systemic function. Since other systemic psfs may be affected 
by such a change - this can be detected by re-computing the weighted dynamic graph - more 
adjustments of elementary rate parameters may be indicated, possibly by an iterative process 
(cf. [33] ). 

Sampling for multiple inputs yields a transfer function matrix, which can be analyzed 
for dependence of the target concentration on each input separately. This can be done by 
standard matrix analysis such as principal component analysis (PCA). For our example, we 
show how species which are poised to integrate signals from two different sources do this under 
the numeric conditions (Fig. [8]). For cAMP production (AC5), we find that AC5GoaGTP 
is dependent only on Da (Fig. [HJ 1), AC5Ca is only dependent on Ca (Fig. [HJ 2), while 
AC5CaGoaGTP is almost not activated at all (Fig. |HJ 3). Even though a link of reactions 
(Ca-AC5Ca-AC5CaGoaGTP-cAMP) exists, signal integration of Ca and Da on AC5 fails 
because of the weak transmission from GoaGTP to AC5Ca. Signal integration between Ca 
and Da occurs for cAMP degradation by calcium-dependent calmodulin regulation of PDE1. 
PP2A with the two variants (calcium-activated) PP2Ac and (PKAc-activated) PP2Ap is 
another potential source of signal integration (Fig. |HJ 4-6) . The psf analysis shows when signal 
integration occurs (here: Da having influence on PP2Ac), and when this effect is negligible 
(here: Ca not having influence on PP2Ap). This may now be studied for correspondence 
with the biological situation. These results emphasize the necessity for numeric analysis of 
input-dependence, beyond the mere existence of links. 

Systemic Delay and State-change Dynamics 

We would like to be able to use systemic psfs with their simple and transparent mathematical 
structure for dynamical simulations. This allows direct experimental testing and fitting by 



10 



bpecies 


U.Ub — > U.bfiM 


U.Ub — > 4.5/iM 


U.o — > U.Ub/UM 


4.5 — > U.Ub/UM 


UaUlrv 


O.D 


U.o 


in 1 

1U.1 


1U 


AL-5*orOalor 1 r 


T C 
< .O 




1 


on 


AbobaLroaLr 1 r 




2.4 


19.5 


26 


cAMP 


7.8 


3.1 


18.1 


18 


PKAc 


251 


164 


347 


300 


pThr34 


288 


200 


369 


307 


PP2Ap 


511 


387 


706 


683 


PP2Ac 


600 


483 


756 


717 


pThr75 


493 


359 


756 


755 



Table 5. For near instantaneous input signal Da at shown concentrations, the table shows 
delays (in s) to reach EC5/EC95 for target species. Decrease is usually slower than 
increase, due to the asymmetry of kon/koff binding parameters. Delay times are sensitive 
both to the amount of increase and the absolute value of the increase, with delays being 
faster for larger signals. 

time series measurements beyond dose-response relationships. In order to do this, we need to 
compute the systemic delays, i.e. the reaction time until a steady state is formed. Then we 
can build a state-change dynamical model from systemic psfs alone, using the appropriate 
delays for the input and the update of a system state. 

Systemic delays depend on the absolute size of the signal and also the direction (in- 
crease/decrease) of signaling. Delays for species in the example system in response to input 
are shown in Table For the computation of target concentrations, we only need a ratio 
such as -j^jj = Kd (binding coefficient) or for forward and backward enzymatic reac- 
tions. For the delays, the difference between kon and koff or kcat and kcatl defines reaction 
times for synthesis and degradation. Therefore, delay computations are fairly complex, but 
the results are often within a fairly narrow range for each reaction (Table [5]). For discrete 
state-change simulations we may use maximal delays for each species. 

From a biological perspective, this table provides an important test on the validity of 
the model. In many cases, systemic delays can be measured. For instance, the delay for 
PKAc at 150-250s rather than 30-60s, as measured in [38] (cf. [7]), seems large and may be 
an indication for a revision of the underlying parameters. From the theoretical perspective, 
this system seems to operate on separate time scales: l-10s, 150-300s and 450-600s. Such a 
separation of reactions by their characteristic delay times is interesting, since it could lead 
to simulation models with different discrete time scales. Here we may calculate psf values 
for fast species with 10s time resolution, for intermediate species with 300s time resolution, 
and for slow species with 600s time resolution, i.e. system update time for state changes 
(Fig. [Til]) . It is an empirical question, whether separate time scales rather than a continuum 
of delay values will prove to be an organizing principle in protein signaling systems [30J. 
A general study, for instance, using models from the BioModels Database (23], might give 
answers to this question. Time scale separation may provide a conservative property of a 
signaling system against fluctuations of concentrations. If total concentrations in the system 
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change, e.g. by protein expression up- or down regulation, miRNA interaction, or diffusional 
processes across compartments, the relevant interactions will continue within each time slice. 
Concerted regulation of protein expression levels may set a clock for the rapidity of signal 
transduction. 

Systemic dynamics, in contrast to elementary reaction dynamics, need not follow a hy- 
perbolic curve. If there are feedbacks in the system, the dynamics may contain transients, 
i.e. the concentration may be higher or lower before it settles into its steady-state value. 
The dynamic response of target species to input are shown in Fig. [9j For a species with- 
out a transient response, the actual value of a species at a shorter delay is always bounded 
by the steady-state value, and all possible concentrations in a continuous-time dynamical 
system are bounded by the psf concentration range [30]. However, if there are transients, 
a psf-based dynamical simulation will miss these transients and plot a simplified trajectory. 
This means that results from a psf analysis with slow inputs cannot be extrapolated to much 
faster input dynamics - in contrast to continuous-timed dynamical systems where arbitrary 
time units can be chosen. This restriction may capture a biological reality: steady-state 
behavior provides the framework and may operate according to rules and principles which 
are separate from the effects of short term fluctuations. 

The psf system allows to generate a dynamical system as a sequence of states defined 
by fluctuations of input. Fig. [TU1 shows an overlay of a differential equation simulation and 
psf state change simulation for a sequence of inputs with 10s duration. Accordingly, the psf 
approximation is excellent for all species with a delay time of < 10s. If we plot psf values for 
slow species at intervals corresponding to their maximal delays, we may linearly interpolate 
between points, and in this case achieve a good approximation of the continuous-time model. 

A psf-based dynamical system is an important tool in order to generate a time-series 
simulation from a calculated model system for comparison with experimental data. The psf 
system utilizes parameters which are uniform and have linear error ranges (cf. Fig.[H Fig. |3j, 
and therefore should improve interaction of the model with the experimental reality. The 
psf model will also allow to predict the optimal stimulation times for different inputs such 
that responses can be measured in steady-state. 

Computing Input-dependent Modularity 

Since we are able to define signal transmission capacity, we have a tool to investigate modu- 
larity of a signaling system. As species saturate or return to basal levels, they act as inactive 
links, i.e. they are stuck at the same concentration value, and cannot transmit further in- 
creases or decreases of input signals. We hypothesize that this effect is actually important 
and dramatic in many protein signaling systems. We may define an inactive connection 
as a species node which has only limited (e.g. < 10%) signal transmission capacity. The 
interconnectedness of a system is then proportional to the number of inactive connections, 
and a module is a part of the system with few or no active connections to the rest of the 
system. 

In the following we discuss the activation/inactivation of links with respect to input 
increases. I.e. given a certain level of extracellular signaling, what happens if this level is 
raised and then kept at the higher level for some time, sufficient for the system to settle into 
a new steady state? Which species nodes will respond to the increase and transmit it to 
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downstream targets, and which species nodes will become saturated and only respond with 
their saturated value? It is also clear that species which have become saturated (inactive) 
will not respond to fast extracellular signals anymore. This analysis is therefore useful both 
for the steady-state context and for understanding of fast input fluctuations. 

For each target species, there is a saturating input level. If that input level has been 
reached, the species can be considered to have become an inactive connection, i.e. a node 
which does not transmit signals. We may define systemic psfs for 'input-target psfs' from the 
input to any target species (Fig. [TT]) . We notice how number of steps in the computation of 
a species concentration relates to a lower cut-off value for signal transmission (e.g. DaDIR, 
AC5GoaGTP, cAMP vs. PKAc, pThr34, pThr75). In other words, earlier steps in the 
sequence saturate at a higher input level than later steps. Two parallel targets (pThr34, 
PP2Ap) of an intermediate step (PKAc) may saturate at very different input levels (~ 3/iM 
for pThr34 vs. 1.5//M for PP2Ap). This mechanism demonstrates the effect of a sequence 
of saturating functions, and constitutes a general principle in the construction of a signal 
transduction system. 

In Fig. [131 we show modularization of the system under Da input with a fixed, fairly high 
concentration of Ca (8/iM). Species which are proximal to Da input, the receptor-ligand 
complex DaDIR, the signaling complex through G proteins and adenylyl cyclase AC5, as well 
as cAMP, are most responsive to Da over a large range of input. Species in the 'integration 
zone' between Da and Ca, such as DARPP-32, PP2A cease responding to Da increases 
at lower levels and become inactive links at higher levels. Species in the system with no 
significant change in concentration at any input level are for the most part proximal to Ca 
instead of Da, and thus highlight modularity among pathways. We are analyzing a biological 
system with two 'pathways', the cAMP and the calcium pathway, which are cross-linked in a 
convergence zone of species which are influenced by both pathways. In this case, we see that 
Da inputs are transmitted to distant targets only up to an intermediate range and that there 
are a significant number of species which do not react to Da at all. Above that range, even 
though closely coupled targets still respond to the input, the signal increase is not registered 
beyond cAMP production and synthesis. 

Similar observations apply when Ca is input (Fig. [T2|) and Da is set to a high value 
(> 1/iM). In this case, the calcium-responsive proteins like calmodulin, CaMKII, calcium- 
activated PP2B (calcineurin) transmit signals, while the GPCR pathway remains almost 
completely unresponsive. There is some signal integration with calmodulin-activated PDE1 
for cAMP. Other than that, we see that PP2Ac and the pThr34 variant of DARPP-32 
respond strongest to calcium while PP2Ap and the pThr75 variant of DARPP-32 have less 
or no responsiveness to calcium. It is remarkable in this case how clearly three different 
modules appear: the GPCR/cAMP pathway, the Ca pathway, and the signal integration 
zone between them. We have kept the second input at a fixed, high value, to simplify the 
graphic for a single variable input. If we combine variable Da with low Ca or variable Ca 
with low Da, (which defines the whole input space), the overall responsiveness of the system 
increases (supplement figures). Because of the specific structure of this system, the two 
input pathways remain segregated, and responsiveness increases specifically in the 'signal 
integration zone' (DARPP-32, PP2A and its variants). 

Signal integration from Ca and Da input is strongest when both Ca and Da are low. This 
is a direct consequence of the effect of coupling saturating nodes. It shows a widespread 
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interconnectedness at low input levels and a stronger modularization with many inactive 
links at higher input levels which results from uniformly saturating activation curves for all 
relevant biochemical reactions. In general, then, modularization, the separation of signaling 
complexes by inactive links, is more pronounced at higher input levels. 

Discussion 

Continuous dynamical systems - systems that use change over time as a system primitive 
- are notoriously difficult to analyze and may not be the best choice of a tool for signal 
transduction systems of moderate or large size. Steady-state matrix computations are sim- 
ple and fast, scale well to very large sizes, and offer multiple opportunities for analysis. By 
calculating transfer functions from a systemic dynamical model, we also gain the opportu- 
nity to extract and analyze parts of the system. This may help in creating re-usable system 
parts. This paper demonstrated a transformation of a mass-action kinetic biochemical re- 
action model implemented by a set of differential equations into an input-response transfer 
function model. The transformation is done by calculating steady-state concentrations for 
each species in response to a range of input values, and then analyzing the resulting vectors 
(matrices) as the basis of a transfer function ('psf'). It could be shown that psfs are also suf- 
ficient to create discrete time dynamical models, with certain restrictions on fast dynamical 
inputs below the time resolution of the system. The primary focus of the analysis was on 
single-input systems, where signaling functions can be matched by parameterized hyperbolic 
functions. However, the analysis can be extended to multiple-input systems by computing 
and analyzing psf transfer matrices. The psf functions allow to cut through the complexity 
of the model and examine interactions in a localized way. If continuous dynamical models 
are being used with the goal of adequately simulating cellular processes, this kind of analysis 
is probably indispensable. The analytical tools of the psf approach may be used to system- 
atically investigate the effect of localized changes to the system [321 133] and to offer local 
corrections to the complex system in a transparent way. 

The extraction of input-target psfs with characteristic hyperbolic saturating properties 
allows to determine input level dependent inactivation of a species node. Accordingly, we 
can define the limits of signal transmission by the distribution of inactive nodes. In the 
example system, a strong distinction of calcium- and cAMP-dependent pathways and a 
signal integration zone were revealed. We could see that at low levels of input, widespread 
interactions are possible, while at higher levels of input, many species enter into a state of a 
constant function value and become inactive links. 

Signal transduction may also be analyzed from the perspective of rational system design. 
Such work is still in its infancy. We may for instance investigate the effect of negative 
feedback links on concentration ranges and times to steady state. Another question would 
be the optimization of a signal transduction system for the trade-off between speed and 
efficacy of signal transmission. With this choice of model, many new questions can be 
raised, and old problems like parameter dependency, modularity or signal integration can be 
addressed in a novel way. 



14 



References 

1. Lindskog M, Kim M, Wikstrom MA, Blackwell KT, Kotaleski JH (2006) Transient 
calcium and dopamine increase PKA activity and DARPP-32 phosphorylation. PLoS 
Computational Biology 2: ell9. 

2. Svenningsson P, Nishi A, Fisone G, Girault J, Nairn A, et al. (2004) Darpp-32: An 
integrator of neurotransmission. Annual Review of Pharmacology and Toxicology 44: 
269-296. 

3. Cepeda C, Levine MS (2006) Where do you think you are going? The NMDA-D1 
receptor trap. Science's STKE : signal transduction knowledge environment 2006: 
pe20. 

4. Yger M, Girault JA (2011) DARPP-32, Jack of All Trades. . . Master of Which? Fron- 
tiers in Behavioral Neuroscience 5: 56. 

5. Kotaleski JH, Blackwell KT (2010) Modelling the molecular mechanisms of synaptic 
plasticity using systems biology approaches. Nature Reviews Neuroscience 11: 239-51. 

6. Stefan MI, Edelstein SJ, Le Novere N (2008) An allosteric model of calmodulin explains 
differential activation of PP2B and CaMKII. Proceedings of the National Academy 
of Sciences of the United States of America 105: 10768-73. 

7. Kim M, Huang T, Abel T, Blackwell KT (2010) Temporal sensitivity of protein kinase 
a activation in late-phase long term potentiation. PLoS Computational Biology 6(2): 
el000691. 

8. Pepke S, Kinzer-Ursem T, Mihalas S, Kennedy MB (2010) A dynamic model of inter- 
actions of Ca2+, calmodulin, and catalytic subunits of Ca2+/calmodulin-dependent 
protein kinase II. PLoS Computational Biology 6: el000675. 

9. Manninen T, Hituri K, Kotaleski JH, Blackwell KT, Linne ML (2010) Postsynaptic 
signal transduction models for long-term potentiation and depression. Frontiers in 
computational neuroscience 4: 152. 

10. Nakano T, Doi T, Yoshimoto J, Doya K (2010) A kinetic model of dopamine- 
and calcium- dependent striatal synaptic plasticity. PLoS Computational Biology 6: 
el000670. 

11. Fernandez E, Schiappa R, Girault JA, Le Novere N (2005) DARPP-32 is a robust 
integrator of dopamine and glutamate signals. PLoS Computational Biology 2(12): 
el76. 

12. Le Novere N, Li L, Girault JA (2008) DARPP-32: molecular integration of phospho- 
rylation potential. Cellular and Molecular Life Sciences : CMLS 65: 2125-7. 

13. Hughey JJ, Lee TK, Covert MW (2010) Computational Modeling of Mammalian Sig- 
naling Networks. Wiley Interdiscip Rev Syst Biol Med 2(2): 194-209. 



15 



14. Qi Z, Miller GW, Voit EO (2008) Computational systems analysis of dopamine 
metabolism. PloS One 3: e2444. 

15. Qi Z, Miller GW, Voit EO (2010) The internal state of medium spiny neurons varies 
in response to different input signals. BMC Systems Biology 4: 26. 

16. Oliveira RF, Kim M, Blackwell KT (2012) Subcellular Location of PKA Controls 
Striatal Plasticity: Stochastic Simulations in Spiny Dendrites. PLoS Computational 
Biology 8: el002383. 

17. Barbano PE, Spivak M, Flajolet M, Nairn AC, Greengard P, et al. (2007) A math- 
ematical tool for exploring the dynamics of biological networks. Proceedings of the 
National Academy of Sciences of the United States of America 104: 19169-74. 

18. Bhalla US, Iyengar R (1999) Emergent Properties of Networks of Biological Signaling 
Pathways. Science 283: 381-387. 

19. Gilbert D, Fuss H, Gu X, Orton R, Robinson S, et al. (2006) Computational method- 
ologies for modelling, analysis and simulation of signalling networks. Briefings in 
bioinformatics 7: 339-53. 

20. Wu J, Vidakovic B, Voit EO (2011) Constructing stochastic models from deterministic 
process equations by propensity adjustment. BMC Systems Biology 5: 187. 

21. Erhard F, Friedel CC, Zimmer R (2008) FERN - a Java framework for stochastic 
simulation and evaluation of reaction networks. BMC Bioinformatics 9: 356. 

22. Bandara S, Schloder J, Eils R, Bock H, Meyer T (2009) Optimal experimental design 
for parameter estimation of a cell signaling model. PLoS Computational Biology 5: 
el000558. 

23. Li C, Donizelli M, Rodriguez N, Dharuri H, Endler L, et al. (2010) BioModels 
Database: An enhanced, curated and annotated resource for published quantitative 
kinetic models. BMC Systems Biology 4: 92. 

24. Deng J, Feinberg M, Jones C, Nachman A (2011) On the steady states of weakly 
reversible chemical reaction networks. Arxiv preprint arXiv:11112386v2 : 1-20. 

25. Lagarias J, Reeds J, Wright M, Wright P (1998) Convergence properties of the nelder- 
mead simplex method in low dimensions. SIAM Journal of Optimization 9: 112-147. 

26. Kraeutler MJ, Soltis AR, Saucerman JJ (2010) Modeling cardiac /3-adrenergic signal- 
ing with normalized-Hill differential equations: comparison with a biochemical model. 
BMC Systems Biology 4: 157. 

27. Yu RC, Resnekov O, Abola AP, Andrews SS, Kirsten R, et al. (2010) The Alpha 
Project: a model system for systems biology research. IET systems biology 2: 222- 
233. 



16 



28. Gutenkunst RN, Waterfall JJ, Casey FP, Brown KS, Myers CR, et al. (2007) Univer- 
sally sloppy parameter sensitivities in systems biology models. PLoS computational 
biology 3: 1871-78. 

29. Theis F, Bohl S, Klingmueller U (2011) Theoretical analysis of time-to-peak responses 
inbiological reaction networks. Bulletin of Mathematical Biology 73: 978-1003. 

30. Purvis JE, Radhakrishnan R, Diamond SL (2009) Steady-state kinetic modeling con- 
strains cellular resting states and dynamic behavior. PLoS computational biology 5: 
el000298. 

31. Snyder GL, Galdi S, Fienberg AA, Allen P, Nairn AC, et al. (2003) Regulation of 
AMPA receptor dephosphorylation by glutamate receptor agonists. Neuropharmacol- 
ogy 45: 703-713. 

32. Maslov S, Ispolatov I (2007) Propagation of large concentration changes in reversible 
protein-binding networks. Proceedings of the National Academy of Sciences of the 
United States of America 104: 13655-60. 

33. Maslov S, Sneppen K, Ispolatov I (2007) Spreading out of perturbations in reversible 
reaction networks. New journal of physics 9: 273. 

34. Famili I, Mahadevan R, Palsson BO (2005) k-Cone analysis: determining all candidate 
values for kinetic parameters on a network scale. Biophysical journal 88: 1616-25. 

35. Zamora-Sillero E, Hafner M, Ibig A, Stelling J, Wagner A (2011) Efficient charac- 
terization of high- dimensional parameter spaces for systems biology. BMC Systems 
Biology 5: 142. 

36. Pozo C, Marin-Sanguino A, Alves R, Guillen-Gosalbez G, Jimenez L, et al. (2011) 
Steady-state global optimization of metabolic non-linear dynamic models through re- 
casting into power-law canonical models. BMC Systems Biology 5: 137. 

37. Smallbone K, Simeonidis E, Swainston N, Mendes P (2010) Towards a genome-scale 
kinetic model of cellular metabolism. BMC Systems Biology 4: 6. 

38. Woo NH, Duffy SN, Abel T, Nguyen PV (2003) Temporal spacing of synaptic stimu- 
lation critically modulates the dependence of LTP on cyclic AMP-dependent protein 
kinase. Hippocampus 13(2):293-300. 



17 




Figure 1. Derivation of a psf for complex formation: PSFs were generated for 
[A] + [B] <B- [AB] with Kd = 200, 500, 2000, [A] = lOnM - 10/iM; [B] = 7fiM, and fitted 
with the saturating hyperbolic function shown in the figure. For an elementary reaction, 
C = Kd. All other curves were generated only by varying C, where ymax=[S]o and n=l. 
The slope at C indicates the signal transmission strength. 




Figure 2. Dynamics of complex formation depends on concentration: 

[A] + [B] <B- [AB] with Kd = 200, [B] = 200nM; [A] = 200nM - 7fiM. Increasing 
concentrations for the source species [A] speeds up the reaction. 
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Figure 3. Enzymatic Reactions with Reverse Reaction. Example reaction with 
[A] = 200, kon=0.00026;koff=1.5; kcat=30.4; kcat2=0.26 (shown in red) and variations 
(shown in green, blue). Both the forward reaction (kcat) and reverse reaction (kcat2) 
parameters are varied by 30%-200%. Variability in signal transmission is expressed by the 
fitted parameter C (half-maximal activation) in a uniform way. 
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Figure 4. Signal Transmission Strength: A 100% increase of input yields diminishing 
increases at higher concentrations. Signals are shown at EC5, EC50 and EC95. 




Figure 5. Weighted Dynamic Network View of the Model System. Input conditions are Ca=8/iM and 
Da=60nM — 5/iM . Each species node is labeled with its concentration range. Reaction nodes are labeled by parameters for a 
hyperbolic fit, or a linear fit if that was sufficient. Complex formation reactions (red) are labeled for both [A] — > [AB] and 
[B] — > [AB] (left to right), enzymatic reactions (blue) are only labeled for [E] —> [A*]. 
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Figure 6. Ranges of Concentrations in Response to Input (Da). For a number of 
relevant target species from the biological model, the range of concentrations for input 
(Da) from lOOnM to 5 fiM is shown (Ca - high, 8 fiM). 
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Figure 7. Local adjustment of a transfer function. A. Elementary (dashed) and 
systemic (continuous) psfs for two targets of GoaGTP. A new systemic psf for 
GoaGTP-AC5CaGoaGTP is defined by functional parameter adjustment (black line). B. 
Elementary parameters were changed to match the new systemic psf. For AC5CaGoaGTP, 
kon=0.0192,koff=25 was adjusted to kon=0.022,koff=1.5, for AC5GoaGTP, 
kon=0.0385,koff=50 was adjusted to kon=0.0495,koff=48.5. 
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D E F 

Figure 8. Dependence of target species on Ca and Da inputs. The 2D psf shows 
species which segregate to only one of the input pathways (A,B,F), a species which remains 
unresponsive (C), and species which show some signal integration in their input 
constellation (D,E). 
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Figure 9. Dynamics of target species to Da input. A. Fast species (<10s) have 
transients. The psf approximation only calculates the steady-state value. B. For slow 
species no transients are apparent. It may take several minutes to reach steady state. 
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Figure 10. Dynamic Response to Da Input. Dots mark systemic psf values, thick 
lines are continuous dynamical simulations by differential equations, thin lines are 
interpolations between 10s psf values. Red arrows mark time points for interpolations for 
slow species like PP2Ac, PP2Ap, pThr34. 
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Figure 11. Systemic psfs for input to target. Shown are EC90 and EC10 
concentrations as cut-off thresholds for effective signaling by input and the slope values at 
threshold. 
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Figure 12. Modularity in signal transmission. Species nodes are colored according to their saturation/depeletion status 
(EC90/EC10) in response to percentage of input. Input is Ca (100nM...10/iM), Da is set at 5 \iM. The figure shows a large 
number of species which have no response (less than 10%) to Ca input, two separate pathways are apparent. Species 
inactivate at low input levels in the 'integration zone'. 




Figure 13. Modularity in signal transmission. Species nodes are colored according to their saturation/depeletion status 
(EC90/EC10) in response to percentage of input. Input is Da (20nM...5/iM), Ca is set at 8 \iM. The graph shows 
input-dependent inactivation of links. Individual effects can be studied. For instance, PDE4 (in contrast to PDEl,PDElCaM) 
shows responsivity to Da input because of a larger complex formation with cAMP, which subtracts from the enzyme 
concentration. 




Figure 14. Graph representation of signal transmission. Input concentration for Da is 20nM . . . 5/jlM, for Ca is 
lOOnM. There is widespread interconnectedness for low inputs. 
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Figure 15. Graph representation of signal transmission. Input concentration for Da is lOOnM, for Ca is lOOnM 
. . . 10/zM. The Da/GPCR pathway is clearly separated from Ca inputs. 



